/* START 3_genHomeWork.do */

/* this script generates home and work locations for the urban/spatial model. We
adopt a week by week approach:
- home: modal origin in the morning/destination in the evening
- work: modal destination in the morning/origin in the evening
- leisure: any weekend trip away from home, divided by the total number */

/* prerequisite packages: fillmissing, labutil2, gtools, parallel */
version 16.1
set more off
cap clear

/* PATHS */
qui do "code/config.do"

/* 20231002: parallel no longer needs this */
/* parallel initialize 8 */

/* load accessory functions */
qui do "./code/99_accessoryFunctions.do"

local panel_flag 0 /* operate only on panel data */

if `panel_flag' == 1 {
	use "make_data/panel.dta", clear
	}
else {
	use "make_data/2015.dta", clear
	cap drop pairs
	}

/* fix labeling from earlier */
cap label copy oa1 oa
cap label drop oa1
cap label values oa da oa

/***************************/
/** SUBSET THE LONG PANEL **/
/***************************/

local planning_areas op dp
local subzones oa da
local units_of_analysis `subzones'
/* local units_of_analysis `subzones' `planning_areas' */
// local units_of_analysis `planning_areas'
local other_restrictions & (TKT_TYP_CD == 34 | TKT_TYP_CD == 38)
local suffix AdultsOnly

/* loop over all trip dates and compute home/work */
gen tripDate = dofc(START)
format tripDate %td

foreach yr in 2015 2016 2017 2018 {
	foreach mth in mar jun sep dec {
		local qr `mth'`yr'
		disp "Now computing `units_of_analysis' ///
		  flows for `qr' with suffix `suffix'"
		cap fd t`qr'
		/* frame put in 1/2, into(t`qr') */
		frame put if (tripDate >= td("01`qr'") & ///
		  tripDate <= td("09`qr'")) `other_restrictions', into(t`qr')
		fv t`qr'
		if (strpos("`units_of_analysis'", "op") > 0) drop if op == . | dp == .
		if (strpos("`units_of_analysis'", "oa") > 0) drop if oa == . | da == .
		if _N > 0 {
			/* execute modal OD code */
			cap gen byte hstart = hh(START)
			cap gen byte am = hstart >= 5 & hstart <= 11
			cap gen byte pm = hstart >= 15 & hstart <= 23
			cap gen dayOfWeek = dow(dofc(START))

			/* figure out correct locations */
			cap restore, not
			cap fd w`qr'
			frame put if dayOfWeek >= 1 & dayOfWeek <= 5, into(w`qr')
			frame w`qr' {
				foreach timeOfDay in am pm {
					/* GENERATE MODAL HOME/WORK */
					foreach x in `units_of_analysis' {
						disp "Now generating `x'_mode_`timeOfDay'"
						cap lab drop `x'_mode_`timeOfDay'
						/* drop if `timeOfDay' == 0 */
						hashsort CRD_NUM `timeOfDay' `x'
						order CRD_NUM `timeOfDay' `x'
						parallel, by(CRD_NUM `timeOfDay'): egen int `x'_mode_`timeOfDay' = ///
						  mode(`x') if `timeOfDay' == 1, ///
						  by(CRD_NUM `timeOfDay') minmode missing
						parallel, by(CRD_NUM): fillmissing `x'_mode_`timeOfDay', by(CRD_NUM)
						}

					/* LABEL VALUES */
					foreach x in `planning_areas' {
						cap label values `x'_mode_`timeOfDay' op
						}

					foreach x in `subzones' {
						cap label values `x'_mode_`timeOfDay' oa
						}
					}

				/************************/
				/* assign home and work */
				/************************/

				foreach name in pHome aHome pWork aWork {
					cap label drop `name'
					}

				/* planning areas only */
				if strpos("`units_of_analysis'", "op") > 0 {
					/* similarly specify home and work for planning areas */
					gen pHome = op_mode_am if op_mode_am == dp_mode_pm
					replace pHome = op_mode_am if dp_mode_am == .
					replace pHome = dp_mode_pm if op_mode_am == . | pHome == .
					replace pHome = op_mode_am if pHome == .
					lab val pHome op
					gen pWork = dp_mode_am if dp_mode_am == op_mode_pm
					replace pWork = dp_mode_am if op_mode_pm == .
					replace pWork = op_mode_pm if dp_mode_am == . | pWork == .
					replace pWork = dp_mode_am if pWork == .
					lab val pWork op

					cap fd work`qr'
					frame put CRD_NUM TKT_TYP_CD pHome pWork, into(work`qr')
					frame work`qr' {
						gduplicates drop
						drop if pHome == . | pWork == .
						gcollapse (count) totwork = CRD_NUM, by(pHome pWork) merge
						cap gcollapse (count) totworkw = CRD_NUM if TKT_TYP_CD == 38, ///
						  by(pHome pWork) merge
						gcollapse (mean) totwork (mean) totworkw, by(pHome pWork)
						ren pWork dp
						}
					}

				/* subzones only */
				if strpos("`units_of_analysis'", "oa") > 0 {
					gen aHome = oa_mode_am if oa_mode_am == da_mode_pm
					replace aHome = oa_mode_am if da_mode_am == .
					replace aHome = da_mode_pm if oa_mode_am == . | aHome == .
					replace aHome = oa_mode_am if aHome == .
					lab val aHome oa
					gen aWork = da_mode_am if da_mode_am == oa_mode_pm
					replace aWork = da_mode_am if oa_mode_pm == .
					replace aWork = oa_mode_pm if da_mode_am == . | aWork == .
					replace aWork = da_mode_am if aWork == .
					lab val aWork oa

					cap fd work`qr'
					frame put CRD_NUM TKT_TYP_CD aHome aWork, into(work`qr')
					frame work`qr' {
						gduplicates drop
						drop if aHome == . | aWork == .
						gcollapse (count) totwork = CRD_NUM, by(aHome aWork) merge
						cap gcollapse (count) totworkw = CRD_NUM if TKT_TYP_CD == 38, ///
						  by(aHome aWork) merge
						gcollapse (mean) totwork (mean) totworkw, by(aHome aWork)
						ren aWork da
						}
				}

				/* NOTE: STILL DOESN'T WORK WITH OP AND OA TOGETHER */
				if (strpos("`units_of_analysis'", "op") > 0) keep CRD_NUM TKT_TYP_CD pHome pWork
				if (strpos("`units_of_analysis'", "oa") > 0) keep CRD_NUM TKT_TYP_CD aHome aWork
				gduplicates drop
				gduplicates drop CRD_NUM, force
				}
			/* fv t`qr' */
			frlink m:1 CRD_NUM, frame(w`qr')
			if (strpos("`units_of_analysis'", "op") > 0) frget pHome pWork, from(w`qr')
			if (strpos("`units_of_analysis'", "oa") > 0) frget aHome aWork, from(w`qr')

			/* preserve */
			/* keep weekend trips only */
			if (strpos("`units_of_analysis'", "op") > 0) {
				keep if (pHome != . & pWork != .) & mod(dayOfWeek, 6) == 0
				gen leisure = pHome != dp | op == dp
				gcollapse (count) dayOfWeek (mean) TKT_TYP_CD, by(CRD_NUM pHome dp)
				hashsort CRD_NUM dp dayOfWeek
				}
			else if (strpos("`units_of_analysis'", "oa") > 0) {
				keep if (aHome != . & aWork != .) & mod(dayOfWeek, 6) == 0
				gen leisure = aHome != da | oa == da
				gcollapse (count) dayOfWeek (mean) TKT_TYP_CD, by(CRD_NUM aHome da)
				hashsort CRD_NUM da dayOfWeek
				}

			ren dayOfWeek leisureTrips

			/* Generate modal leisure destination */
			gegen maxTrips = max(leisureTrips), by(CRD_NUM)
			gen argmax = 1 if leisureTrips == maxTrips
			/* gegen numArgmax = total(argmax), by(CRD_NUM) */

			frames put CRD_NUM TKT_TYP_CD *Home d* argmax, into(mode`qr')
			frame mode`qr'{
				keep if argmax == 1
				hashsort CRD_NUM
				by CRD_NUM: sample 1, count
				/* if repeated max, draw at random. DESTRUCTIVE STEP */
				if (strpos("`units_of_analysis'", "op") > 0) {
					gcollapse (count) totleis = CRD_NUM, by(pHome dp) merge
					cap gcollapse (count) totleisw = CRD_NUM if TKT_TYP_CD == 38, ///
					  by(pHome dp) merge
					gcollapse (mean) totleis (mean) totleisw, by(pHome dp)
					keep if dp != .
					/* fill missing planning areas */
					xtset pHome dp
					tsfill
					xtset dp pHome
					tsfill
					}

				else if (strpos("`units_of_analysis'", "oa") > 0) {
					gcollapse (count) totleis = CRD_NUM, by(aHome da) merge
					cap gcollapse (count) totleisw = CRD_NUM if TKT_TYP_CD == 38, ///
					  by(aHome da) merge
					gcollapse (mean) totleis (mean) totleisw, by(aHome da)
					keep if da != .
					/* fill missing subzones */
					xtset aHome da
					tsfill
					xtset da aHome
					tsfill
					}
				}

			/* leisure measure 2: 1/N */
			drop maxTrips argmax
			gegen totLeisCard = total(leisureTrips), by(CRD_NUM)
			gen leisFrac = leisureTrips/totLeisCard

			/* Generate leisure and workfare alt leisure counts */
			if (strpos("`units_of_analysis'", "op") > 0) {
				gcollapse (sum) totleisalt = leisFrac, by(pHome dp) merge
				cap gcollapse (sum) totleiswalt = leisFrac if TKT_TYP_CD == 38, ///
				  by(pHome dp) merge
				cap gcollapse (sum) totleiswcount = leisureTrips if TKT_TYP_CD == 38, ///
				  by(pHome dp) merge
				gcollapse (mean) totleisalt (mean) totleiswalt ///
				  (sum) totleiscount = leisureTrips (mean) totleiswcount, by(pHome dp)
				keep if dp != .
				/* fill missing planning areas */
				xtset pHome dp
				tsfill
				xtset dp pHome
				tsfill
				frlink m:1 pHome dp, frame(work`qr')
				frlink m:1 pHome dp, frame(mode`qr')
				}

			else if (strpos("`units_of_analysis'", "oa") > 0) {
				gcollapse (count) totleisalt = leisFrac, by(aHome da) merge
				cap gcollapse (count) totleiswalt = leisFrac if TKT_TYP_CD == 38, ///
				  by(aHome da) merge
				cap gcollapse (sum) totleiswcount = leisureTrips if TKT_TYP_CD == 38, ///
				  by(aHome da) merge
				gcollapse (mean) totleisalt (mean) totleiswalt ///
				  (sum) totleiscount = leisureTrips (mean) totleiswcount, by(aHome da)
				keep if da != .
				/* fill missing subzones */
				xtset aHome da
				tsfill
				xtset da aHome
				tsfill
				frlink m:1 aHome da, frame(work`qr')
				frlink m:1 aHome da, frame(mode`qr')
				}

			frget totwork totworkw, from(work`qr')
			frget totleis totleisw, from(mode`qr')
			drop work`qr' mode`qr'

			foreach var in totleis totleisw totwork totworkw ///
			  totleisalt totleiswalt totleiscount totleiswcount {
				replace `var' = 0 if `var' == .
				}
			gen int dt = td(01`qr')
			format dt %td
			compress

			/* clean up: delete the non-t frames */
			cap fd work`qr'
			cap fd mode`qr'
			cap fd w`qr'
			if (strpos("`units_of_analysis'", "op") > 0) ///
			  save "make_data/`qr'`suffix'.dta", replace
			if (strpos("`units_of_analysis'", "oa") > 0) ///
			  save "make_data/a`qr'`suffix'.dta", replace
			} /* _N > 0 check ends here */
		frame change default
		}
	}

/* append all to new frame */
cap fd collated
fv collated
foreach yr in 2015 2016 2017 2018 {
	foreach mth in mar jun sep dec {
		local qr `mth'`yr'
		if (strpos("`units_of_analysis'", "op") > 0) {
			cap append using "make_data/`qr'`suffix'.dta"
			/* cap rm "make_data/`qr'.dta" */
			}

		if (strpos("`units_of_analysis'", "oa") > 0) {
			cap append using "make_data/a`qr'`suffix'.dta"
			/* cap rm "make_data/a`qr'.dta" */
			}
		}
	}

if (strpos("`units_of_analysis'", "op") > 0) {
	keep pHome dp tot* dt
	gegen int id = group(pHome dp)
	xtset id dt
	save "make_data/panelODByWeek`suffix'.dta", replace
	}

if (strpos("`units_of_analysis'", "oa") > 0) {
	keep aHome da tot* dt
	cap drop id
	gegen long id = group(aHome da)
	xtset id dt
	save "make_data/panelODAByWeek`suffix'.dta", replace
	}
